Superfluid Insulator Transitions of Hard-Core Bosons on the Checkerboard Lattice 
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We study hard-core bosons on the checkerboard lattice with nearest neighbour unfrustrated hop- 
ping t and 'tetrahedral' plaquette charging energy U. Analytical arguments and Quantum Monte 
Carlo simulations lead us to the conclusion that the system undergoes a zero temperature (T) quan- 
tum phase transition from a superfluid phase at small U/t to a large U/t Mott insulator phase with 
p — 1/4 for a range of values of the chemical potential fx. Further, the quarter-filled insulator breaks 
lattice translation symmetry in a characteristic four-fold ordering pattern, and occupies a lobe of 
finite extent in the fJrU/t phase diagram. A Quantum Monte-Carlo study slightly away from the tip 
of the lobe provides evidence for a direct weakly first-order superfluid-insulator transition away from 
the tip of the lobe. While analytical arguments leads us to conclude that the transition at the tip 
of the lobe belongs to a different landau-forbidden second-order universality class, an extrapolation 
of our numerical results suggests that the size of the first-order jump does not go to zero even at 
the tip of the lobe. 

PACS numbers: 75.10.Jm 05.30.Jp 71.27.+a 



I. INTRODUCTION 

Experiments on a variety of condensed matter systems 
such as under-doped high-temperature superconductors, 
strongly correlated quasi-two dimensional organic com- 
pounds like K-(ET)2Cu2(CN);ji and the triangular lat- 
tice antiferromagnet CsCuCLf^ all point to the possi- 
bility of realizing genuinely new phases of matter out- 
side of standard paradigms such as fermi liquid theory 
of normal metals, spin- wave theory of magnetically or- 
dered states, and the BCS theory of superconductivity. 
There has been a great deal of activity in the recent past 
aimed at developing theoretically consistent descriptions 
of such genuinely new phases of matter — examples in- 
clude various proposals for topologically ordered liquid 
states of strongly-correlated systems, and their gauge- 
theoretic description 3 -^ 5 -! 6 -! 7 -! 8 .! 9 -. 

A generic feature of these 'exotic' phases is that their 
low lying excitations are described in terms of quasipar- 
ticles that carry unusual quantum numbers and are best 
thought of as fractions of the usual electron or magnon 
excitations. These fractionalized quasiparticles interact 
with each other by emergent gauge forces, and have an 
independent existence only in such exotic phases where 
the gauge forces scale to zero at large distances effectively 
'deconfining' the fractionalized quasiparticles. 

In closely related work, Senthil et. al^ have also de- 
veloped a theory of a class of second order quantum phase 
transitions that fall outside of the rubric of Landau's the- 
ory of phase transitions. In the usual Landau theory, 
the transition is described by analyzing fluctuations of 
the order parameter field that encodes the sharp distinc- 
tion between the two phases. In certain situations where 
the phases on either side of the transition break differ- 
ent symmetries, Senthil et. al. demonstrated that this 
usual description fails, in that it would erroneously rule 



out the possibility of a direct second order transition in 
the generic case. 

Such a direct transition is actually possible, and its 
properties are best understood in terms of new frac- 
tionalized degrees of freedom (that can be thought of 
as fractions of the quasiparticle excitations of the ad- 
joining phase) that are liberated at precisely the critical 
ponrfcii. This has sparked interes t 12 ' 13 : 14 ! 15 : 16 in theo- 
retical studies of relatively simple model systems where 
one can check for the presence of such unusual phases 
of matter and non-landau 'deconfined' quantum phase 
transitions with numerical studies. 

In recent worki 7 -, two of the present authors argued 
that a S = 1 antiferromagnet with isotropic exchange 
J and strong uniaxial single-ion anisotropy D on the 
kagome lattice would support an unusual spin-nematic 
state in zero field, and a spin-density wave magneti- 
zation plateau state for a range of non-zero magnetic 
fields along the easy-axis. Both phases and the transi- 
tion between them can be modeled by a Hamiltonian for 
hard-core bosons with unfrustrated hopping t and frus- 
trated nearest-neighbour repulsive interactions U— ; the 
nematic state corresponds to a ordinary superfluid state 
of bosons at half-filling, while the magnetization plateau 
state corresponds to a 1/3 filled charge-density wave in- 
sulator that occupies a lobe in the fi-U ft phase diagram 
(where /i, the chemical potential for bosons, is linearly 
related to the applied magnetic field, and U/t ~ D 2 /J 2 ). 

As the two phases break different symmetries, the tran- 
sition between them has the possibility of being an un- 
usual non-Landau transition 18 . Numerical studies of the 
transition close to the tip of the lob o 17 : 19 reached the 
tentative conclusion that the transition between the two 
phases may indeed be a direct second-order transition 
with unusual exponents, although more recent work on 
larger sizes seems to see faint signatures of a very weak 
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first-order transition close to the tip 1 ^. 

II. THE MODEL AND SCHEMATIC PHASE 
DIAGRAM 

Motivated in part by this intriguing transition seen 
on the Kagome lattice, and its possible relevance to the 
physics of anisotropic antiferromagnets, we consider here 
a closely analogous model of hard-core bosons on the 
checkerboard (planar pyrochlore) lattice, with the aim of 
shedding further light on the physics of such systems (As 
our work was nearing completion, other work on similar 
fermionic model o 20 ' 21 i 22 and somewhat different SU(2) 
symmetric spin model o 23 ' 24 has also appeared in the lit- 
erature). 

The checkerboard or planar pyrochlore lattice is made 
up of four site plaquettes which can be thought of as pro- 
jections of a tetrahedron onto a plane, and are centered 
at the sites of an underlying square net (Fig [1]) . The 
bosonic particles live on the vertices of this 'tetrahedron' 
(which lie on the centers of the links of the square net). 
The Hamiltonian for our boson hubbard model reads 

H=-tJ2(4bj + h.c.) 

<y) 

+ U ^2(nip + n 2p + n 3p + n 4p - f) 2 (1) 

p 

where (ij) refers to the nearest neighbour links of the 
checkerboard lattice, b\ (&,) creates (destroys) a parti- 
cle on site i, t > denotes the nearest neighbour hop- 
ping, and U > denotes the repulsion energy between 
bosons on the same 'tetrahedraP plaquette p. Note that 
n\ p , ri2 P , and n 4p are the boson number operators of 
the sites of the p th such tetrahedral plaquette (which 
can equally well be indexed by the site of the under- 
lying square lattice around which it is centered) , and the 
parameter / serves as a reduced chemical potential that 
tunes the particle density in the ground state. 

Before we begin a detailed discussion of our numeri- 
cal results, it is useful to delineate the gross features of 
the phase diagram using fairly general and reliable argu- 
ments. To this end, we begin in the limit of vanishingly 
small repulsive interaction U -C t. As U/t — » 0, the 
hopping term in the Hamiltonian dominates and its un- 
frustrated nature immediately leads us to the conclusion 
that the ground state is a featureless superfluid. The 
opposite limit when U/t — > oo and the potential energy 
term dominates is best understood by first going over into 
the 'atomic' limit of isolated tetrahedra (t — 0) and then 
analyzing the properties of the perturbation expansion in 
t/U. In the atomic limit, the density on each tetrahedron 
is controlled by the reduced chemical potential /: 

For < / < 1/2, the ground state has no particles, 
while for 1/2 < / < 3/2, the preferred density is one 
particle per tetrahedron (which translates to a density of 




FIG. 1: (Color online), a) Periodic Checkerboard lattice 
and the square lattice whose bonds pass through the checker- 
board sites. The basic corner-sharing units of this lattice are 
'tetrahedral' plaquettes centered at the sites of the square 
net (shown in b)). The presence of a boson on a site of the 
checkerboard lattice can be represented as a 'dimer' cover- 
ing the corresponding link of the underlying square lattice. 
In the plaquette ordered state, red squares resonate via the 
ring-exchange process (shown in c)) while the other links are 
empty (n% = 0). The dotted bonds indicate higher kinetic en- 
ergy in the plaquette state. In the alternate columnar state, 
dimers cover all green edges but not red ones. 



p = 1/4), and for 3/2 < / < 5/2, each tetrahedral pla- 
quette is half-filled. Further, the physics for / > 2 can 
be mapped exactly to the physics at 4 — / by a canon- 
ical transformation (specific to hard-core bosons) which 
is best written in terms of the equivalent pseudospin-1/2 
operators S x - iS y = b, S z = n - 1/2: S y -> -S y , 
S z — » —S z , S x — > S x . Thus the density as a function of 
/ consists of a series of plateaus with a preferred number 
of particles per tetrahedral plaquette, and in view of the 
above, it suffices to consider the physics of the p = 1/4 
and p = 1/2 plateaux to obtain a complete account of 
the physics. 

A crucial feature of these plateaus (except the triv- 
ial cases with no particles and a fully filled lattice) is 
that they represent properties of a manifold of degen- 
erate ground states with macroscopically large degener- 
acy. For concreteness, we focus here on the plateau with 
one particle per tetrahedron (our results for the p = 1/2 
plateau will be discussed separately). A convenient rep- 
resentation of the ground state manifold is to associate 
a dimer on a link of the underlying square lattice with 
the presence of a particle on the corresponding site on 
the checkerboard lattice. Each state in the ground state 
manifold then corresponds to a perfect dimer cover of the 
underlying square lattice. 

The t term introduces quantum tunneling between 
these dimer states, which lifts the classical degeneracy. 
In the dimer language, the leading order effective Hamil- 
tonian in the U/t —* oo limit is a quantum dimer model 
with a ring -exchange term (Fig [J) of 0(t 2 /U) that flips 
two parallel horizontal (vertical) dimers on any flippable 
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t/U 

FIG. 2: Schematic phase diagram for the transition from su- 
perfluid state to Mott insulator state at p = 1/4. Scans I and 
II are discussed in text. 



plaquette to vertical (horizontal) dimers, and acts en- 
tirely within the classical ground state manifold. 

Quantum dimer models with ring-exchange on indi- 
vidual plaquettes of two dimensional bipartite lattices 
(here a square lattice) quite generally have crystalline 
ground states in which the spatial arrangement of dimers 
break lattice symmetry^. In the square lattice case, 
one expects either a plaquette or columnar type Valence 
Bond Solid (VBS), in which there is crystalline order at 
Qi = (it, 0) and Q 2 = (0, n) 26 (lattice directions defined 
in Fig[T]). Thus, at large but finite U/t, one has a Mott 
insulating lobe in the f — U/t phase diagram with den- 
sity pinned to p — 1/4 for a finite range of values around 

As has been discussed earlier by Fisher et. al. in 
their seminal analysis 27 of superfluid insulator transitions 
in bosonic lattice models, the contour of fixed density 
p = 1/4 is expected to terminate on the tip of the lobe, 
while all contours with density different from this value 
skirt around the Mott insulating lobe. Consequently, the 
tip of the lobe of the p = 1/4 Mott insulator represents 
a special point along the locus of superfluid - Mott in- 
sulator transitions — indeed, the low-energy theory at the 
tip is expected to have an additional particle-hole sym- 
metry arising from identical excitation energies for quasi- 
particle and quasi-hole excitations. 

This sort of reasoning leads to the schematic phase di- 
agram shown in Fig [2j but is silent on the nature of the 
quantum phase transition lines between the Mott insulat- 



ing lobes and the superfluid state. Since one of the phases 
is expected to break spatial lattice translation symmetry, 
and the other superfluid state breaks the internal U(l) 
symmetry associated with particle number conservation, 
a Landau-type description of the competition between 
these two very different order parameters would generi- 
cally predict (along the p = 1/4 contour) either a first 
order phase transition, or an intermediate regime with 
both phases or an intermediate regime with neither or- 
der. 

An intermediate phase with neither order would cor- 
respond to a liquid phase of bosons in two dimensions at 
zero temperature. As no consistent scenario for such a 
phase is known in two dimensions in the present context, 
we reach the tentative conclusion that such a phase is 
unlikely. An intermediate phase with both orders would 
correspond to a lattice supersolid, and several examples 
of this are known in the literatur o 28 i 29 i 30 i 31 i 32 i 33 . A direct 
second order phase transition between the two broken 
symmetry phases would require "fine-tuning" to a multi- 
critical point in a Landau type description. However, as 
has been shown recently^ by Senthil et al, conventional 
Landau theory can fail in certain closely related bosonic 
models in which quantum mechanical Berry phase effects 
produce a direct second order phase transition. 



III. NUMERICAL STUDY 

Our goal here is to study this model more closely using 
the well-documented stochastic series expansion (SSE) 
Quantum Monte Carlo (QMC) method 3 ^^ to access 
the phase diagram (Fig [5]) , and then compare our nu- 
merical results with the expectations built on the previ- 
ous work of Senthil et. al. For this purpose, we have 
performed a detailed study along two scans in the U/t- 
f phase diagram: Scan I varies U/t at a fixed value of 
/ — this value has been chosen after several preliminary 
exploratory scans to ensure that it intersects the locus of 
phase transitions close to the tip of the lobe. The sec- 
ond scan (II) was performed at fixed U/t such that one 
encounters an insulating phase for an appreciable range 
of /, and can probe the nature of the ordered state in 
somewhat more detail at a point relatively deep in the in- 
sulator. In addition, to probe nature of the particle-hole 
symmetric transition in more detail, we have also per- 
formed a detailed study of (crystalline) order-parameter 
statistics for a sequence of critical points approaching the 
tip of the lobe. 

Most of our data is on L x L samples (where L x L is the 
number of unit cells) where L ranges from 20 to 48. We 
use standard SSE estimators to calculate the superfluid 
density p s , the equal time correlator of the density S Q/ g(q) 
= (fw(q)^/3r( — q)) an d the static structure factor of the 
density at zero frequency Xa/3(q) = (/ rf , rn QT (q)n / 3o(— q)) 
where n aT (q) = (l/£ 2 )5i i n aT (rj)ea;p(iqr i ), n QT (r 4 ) be- 
ing the boson number operator at imaginary time r, 
and all site types (a — 0,1) are assigned the coordi- 
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FIG. 3: (Color online). Vertical scan (II) showing mott insu- 
lator state at U/t = 3.50 and T = i/10 for size L = 24. Here 
m 2 = \ip c \ 2 , with ip c defined as in Eqn ([2}. Note that the 
magnitude of the bragg peak in the kinetic energy correlator 
is also non-zero in the plateau state; this is exemplified by the 
data shown here (KEoo) on the bragg peak in the correlator 
of kinetic energy along bond-type labeled in Fig [T] Error 
bars are smaller than the symbol sizes in all the cases. 



nates of the corresponding site (i) of the square net 
(Fig [1]). We also calculate the static correlators for 
the 'kinetic energy' K\ = (b\bj + h.c.)i on the link 
I, S^(q) = (JdTK aT (q)Kl (q)) 7 where K aT (q) = 
(V^ 2 ) J2i K aT (ri)exp(iqri), K aT (r) being the kinetic 
energy operator at imaginary time r (all bond types (a — 
0,1,2,3) are assigned the coordinates of the corresponding 
site (i) of the square net as in Fig |T|) . 

To verify the presence of the Mott insulator lobe in the 
phase diagram, we performed vertical scans at constant 
U/t for different values of /. One particular scan (Scan 
II) at U/t = 3.5 is shown here (Fig [3]). The scan clearly 
shows the existence of the mott insulator over an appre- 
ciable region in the phase diagram. As is clear from Fig [4] 
the insulating phase has Bragg peaks in the structure 
factor Soo(q) [^nCci)] anci the static susceptibily xoo(q) 
[Xn(q)] at wavevectors (ir, 0) [(0, w)] which tend to a non- 
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FIG. 4: (Color online). Finite size extrapolations establishing 
the presence of superfluid and charge density wave ordered 
phases. Here / = 0.875 and T = t/15. Error bars are smaller 
than the symbol sizes. 



zero value in the the L — > oo limit. In contrast, there 
are no appreciable Bragg peaks in the cross-correlator 
or off-diagonal static susceptibility Sor (q) [xor(q)]- The 
wavevectors of these peaks point to the presence of either 
plaquette and columnar order in the insulator. 

To obtain a better characterization of this ordering, we 
measure two complex order parameters sensitive to the 
spatial symmetry breaking 



Ipk 



n 0T (Qi) + ini T (Q 2 ) 

{K 2t (Qi) + Kzr{Qi) - K 0t {Q x ) - K lT {Q x )) 
i{K lT (Q 2 ) + K 2T (Q 2 ) - K 0T {Q 2 ) - K 3r {Q 2 )) 



(2) 



where Q\ — (tt, 0) and Q 2 = (0, n). 

This spatial symmetry breaking can now be analyzed 
using a Landau theory written in terms of the order pa- 
rameters ip c and ipk' F° r our purposes here, it suffices 
to focus on just the "potential energy" part of the lan- 
dau theory, and leave out all terms involving spatial gra- 
dients. Terms in the landau theory are constrained by 
the requirement of invariance under the symmetries of 
the system, and the relevant symmetry operations in our 
case are the discrete operations of the space group of the 
square lattice, i.e. tt/2 rotations, translations, reflections 
and inversion. The action of these on the order parame- 
ters is given in Table U 
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Operation 



vr/2 



R 

T 

r, 
n,. 

TZy 
1 



Coordinates 



— » £ijX 



Xi 

y -> 7/ + 1 

a; — > x- + 1 
x — > — x 

y ^ -y 

r — > — r 



i/» c transformation 



V> c — > i^c 

V>c -> ^ 

■<Pc -» 

V>c -> ^ 
V>c — » — tpc 



ipk transformation 



V'fc -> 

V> fc -> -Vfc 
V'fc -> ^ 

— > -Tpk 



TABLE I: Transformations of the order parameters ip c and ^k 
on the square lattice under the discrete symmetry generators 
of the lattice. Here i = 1,2 = x,y is a spatial index. 
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FIG. 5: (Color online). Histograms of the relative phase of 
the two order parameters and the absolute phase of the kinetic 
energy order parameter characterizing the insulator state at 
U/t = 3.50, / = 0.875 and T = t/15 for size L = 28. 



Writing all terms consistent with the symmetries of the 
model (up to fourth order), we obtain the Landau energy 
functional 



Lpot — 



F (|^ c | 2 , |^ fe | 2 ) + CoM + C 4 ) + CeMt + r k 4 ) 

(3) 



where F(x, y) is an analytic function of its arguments. 
The nature of ordering (plaquette versus columnar) is 
determined by the signs of the coefficients C of terms in 
the landau free energy which fix the relative as well as 
absolute phases of the two order parameters. Our results 
on the statistics of the phases of the various order param- 
eters (shown in Fig [5] for a point deep in the insulating 
phase) can be modelled by taking Cg ak to be negative. 
However, the signs of Cg e and Ce k cannot be inferred 
from these histograms as no statistically significant bias 
is apparent in the value of the absolute phase. 

As the absolute phase is too weakly pinned to usefully 
discriminate between the plaquette and columnar order- 
ing possibilities, we need to rely on indirect evidence to 



obtain the nature of the ordering pattern: In this re- 
gard, the presence of clearly discernible Bragg peaks in 
the kinetic energy correlators that survive in the thermo- 
dynamic limit is suggestive of plaquette type order (cor- 
responding to a positive sign for the coefficients Cg c and 
Ce k ) since the simplest caricatures of columnar order- 
ing have essentially no kinetic energy ordering associated 
with the spatial symmetry breaking. Indeed, the nature 
of Bragg peaks seen in the real part of the static corre- 
lators of kinetic energy, S^f{q), is completely consistent 
with expectations based on the simplest variational wave- 
function for the plaquette insulator. More specifically, 
the plaquette ordered variational wavefunction denoted 
pictorially in Fig [1] leads directly to the following be- 
haviour: Re(S™(q)), Re(S%($), Re(S 2 K 2 (q)), Re(Sf(q)) 
have maxima at (tt, 0) and (0, tt), Re(S^{q)), Re{Sj?{q)) 
have maximum at (tt, 0) and minimum at (0,7r), 
Re{S^(q)), Re(Sj((q)) have minimum at (tt, 0) and max- 
imum at (0, tt) and Re(Sj^(q)), Re(S}f(q)) have minima 
at (tt, 0) and (0, tt) (the other correlators can be obtained 
from the above ones trivially). Precisely this behaviour 
is seen in the QMC data shown in Fig [Sj 

In the superfiuid phase, there are no Bragg peaks that 
survive the thermodynamic limit, even in the immediate 
vicinity of the transition to the insulator. However, equal 
time correlator and static structure factor both show 
an interesting 'dipolar' structure in the superfiuid phase 
near the transition (Fig [7]) . These dipolar correlations 
do not seem to be associated with the quantum phase 
transition itself, and persist into the insulating state as 
well. It is therefore instructive to look more closely at the 
static susceptibilities xoo(q) and Xn( c l) an d ask whether 
the dipolar nature is just due the strong local constraint 
introduced by a large U term. 

More specifically, if we are at / ~ 1 and temperatures 
high enough to wash out the coherent quantum dynamics 
associated with the kinetic energy term, one would expect 
static density correlators to match (up to prefactors) the 
predictions for equal time quantities of a classical dimer 
model on the square lattice (viewing each particle as a 
dimer on the corresponding link of the underlying square 
lattice). The dipolar structure seen at lower tempera- 
tures can then be viewed as remnants of this essentially 
classical physics as the system crosses over into a coher- 
ent quantum regime. 

To check this scenario, we fit the static structure factor 
to a functional form closely related to the classical result 
expected for the equal time correlator of classical dimers 
on the square lattice: 



xoo(q) = 
xn(q) = 



a(cos 2 (q y /2) + A 2 ) 
cos 2 (q x /2) + cos 2 (q y /2) + A 2 

q(cos 2 (g x /2) + A 2 ) 
cos 2 (q x /2) + cos 2 (q y /2) + A 2 



(4) 



For the critical state of fully-packed dimers, one expects 
A 2 = 0, and we allow for a non-zero A 2 to account for 
slight density deviations from quarter filling. 
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FIG. 6: (Color online). Plots of Re(S^{q)),Re(Sf{q)), 
Re{S°K(q)) and Re{S%{q)) for L = 28, U/t = 3.190, / = 
0.875 and T = t/15. The presence of the extrema at wavevec- 
tors (n, 0) and (0, it) with the correct signs is consistent with 
the distribution of K.E. shown in Fig[T] Error bars are smaller 
than the symbol sizes. 



O 
O 



0.00065 



0.00045 



0.00025 



5e-05 



o 
o 



0.0006 



0.0004 



0.0002 



(q.«) 
(n.q) 

(q.q) 



▲ 

■ 

T 



-1.5 -1 -0.5 0.5 1 1.5 
q - 71 

L = 28, U/t = 4.0, f = 0.875, T = t 




FIG. 8: (Color online). Top panel: The static susceptibility 
Xoo(q) plotted along three cuts in momentum space for L — 
28, U/t = 3.180, / = 0.875 and T = t/15. Bottom panel: The 
same plotted at U/t = 4, / = 0.875 and higher temperature 
T = t, showing an excellent fit to the dipolar form Eqn (|4]) 
with parameters a — 0.00033 and A 2 = 0.0388. Error bars 
are smaller than the symbol sizes. 
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FIG. 7: (Color online). Static susceptibilities Xoo(q) and 
Xn(q) for L = 28, U/t = 3.180, / = 0.875 and T = t/15. 
In both panels q x and q y vary from to 27r. The white re- 
gions are because values of xoo and xn > 0.001 are set to 
zero. 



As is clear from Fig [51 the dipolar structure in the 
higher temperature data fits remarkably well to our func- 
tional form, with a reasonably small value of A 2 w 0.0388 
reflecting a slight deviation of the density from quar- 
ter filling. Furthermore, consistent with our expectation 
above, the low temperature {fit — 15) data of Fig [5J also 
displays qualitatively similar features, although the clas- 
sical form (EqndJ) does not fit the data as well. 

We now turn our attention to the nature of the quan- 
tum phase transition at the tip of the p = 1/4 lobe. We 
perform Scan I (Fig|2J) at a constant aspect ratio of /3/L 
— 15/28. We also compare data at different sizes at a 
range of inverse temperatures (3 ranging from 5 to 15 to 
ascertain whether the T — * limit is reached for a given 
size. The value of / = 0.875 for this scan was chosen 
after some trial and error to satisfy two criteria: Firstly, 
this appears to be close to the value at which the insu- 
lator has maximum extent. Secondly, the system is close 
to being particle-hole symmetric in the transition region 
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FIG. 9: (Color online). Histograms of the density deviation 
from quarter-filling in the transition region of scan I for size 
L — 28. Deviations are quantified in terms of the total num- 
ber of particles by which the system deviates from quarter- 
filling, and the different curves are labeled by the correspond- 
ing values of U/t at constant / = 0.875 and T — t/15. 



for this value of / (see Fig [SJ) • 

To determine the location of the critical point quantita- 
tively (under the assumption of a scale invariant second- 
order critical point), we calculate the Binder cumulant 
g m and p s L for different system sizes at a constant as- 
pect ratio of ft/L = 15/28. For a continuous phase tran- 
sition, scaling predicts that these two quantities may be 
written as g m = F gm (L 1 / ,/ ((t/U) c - (t/U)), (3/L z ) and 
p s L z = F ps {L l l v Ut/U) c - {t/U)),l3/L z ) near the criti- 
cal point at sufficiently low temperature, where F 9m and 
Fp 3 are the appropriate scaling functions. Hence if the 
transition is continuous with z = 1, both the quantities 
should exhibit the same crossing point in data taken for 
various L at a fixed aspect ratio (3/ L when plotted versus 
U /t. The reasonably sharp crossings and their locations 
indicate that there may be a continuous transition at 
U/t « 3.178 ± 0.002 (to be precise, we obtain a transi- 
tion at U/t w 3.180 based on the crossing of p s L, while 
the binder cumulant of the order parameter m 2 yields a 
crossing point at U/t « 3.176 (see FigQTJ])). The crossings 
shift very little when we consider the constant tempera- 
ture runs at T — t/15 (suggesting that we have reached 
the asymptotic low temperature regime at these sizes). 
However, it appears finite-size effects may be more sig- 
nificant, since we do notice a slight tendency of the super- 
fluid and density wave transitions to move closer to each 
other if we look at the crossings of data at successively 
larger sizes. Our conclusion based on available data at 
these smaller sizes is thus in favour of a single transition 
near the tip of the lobe, with no appreciable region of 
phase coexistence. 

To further explore the possibility that the transi- 
tion is actually continuous and not first order, we do 
a scaling collapse of p s L, x(Q) an d S(Q) (see Fig ITTl). 
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FIG. 10: (Color online). Crossing of the Binder cumulant g m 
= 1 — (m 4 )/(3(m 2 ) 2 ) where m 2 = \4> c \ 2 , and p s L for different 
system sizes at / = 0.875 and j3/L — 15/28. Error bars are 
smaller than the symbol sizes. 



where^x(Q) = (xoo(Qi) + Xu(Q 3 ))/2 and S(Q) = 
(SooiQi) + <5ii(Q2))/2 (the averaging is done to im- 
prove statistics). For a direct continuous transition, 
X (Q) = L-*F x (LV»((t/U) c - (t/U)),(3/L z ) and S(Q) 
= L- 1 i- z F s {L 1 / lJ {{t/U) c - {t/U)),l3/L z ) where the crit- 
ical exponents are the same as in the p s L z scaling rela- 
tion. These scaling relations can also be used to check 
whether z = 1 by extracting 77 from the scaling collapse 
of x(Q) an d then using it for collapsing the S(Q) data. 
From the scaling collapse of the data, we estimate that 
v w 0.46 ± 0.05 and r\ w —0.25 ± 0.12 where the error 
bars are obtained by attempting scaling collapse of the 
data with different values of the exponents. The scaling 
collapses are certainly consistent with a direct transition 
with z = 1, especially since the value of v obtained from 
the superfluid and bragg peak data very nearly coincide 
within their error bars. 

While our scaling collapse does seem reasonable, the 
negative value of r] obtained from these fits leads us to 
speculate that the transition for this particular value of 
/ = 0.875 may actually be very weakly first order — this is 
because a negative value of 77 is ruled out on quite general 
grounds at any conformally invariant critical point. To 
confirm this suspicion, we look at the histograms of the 
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FIG. 11: (Color online). Scaling collapse of the static struc- 
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/ = 0.875 and J3/L — 15/28. Error bars are smaller than the 
symbol sizes. 



charge density wave order parameter to 2 . They show a 
characteristic double-peaked structure when one goes to 
large system sizes indicating that the transition is weakly 
first order (see Fig [T2|). Thus, based on available evi- 
dence, we conclude that the transition close to the tip of 
the lobe is very weakly first order. The first order na- 
ture becomes more apparent at larger sizes L > 48. It 
should be noted that it is precisely at these large sizes 
that the deviation from particle-hole symmetry at crit- 
icality also becomes more evident. This correlation in 
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FIG. 12: (Color online). Distribution of the order parameter 
m 2 for different system sizes at U/t — 3.1775, / = 0.875 and 
P/L = 15/28. The double-peaked structure becomes more 
prominent at larger system sizes, as does the deviation from 
particle-hole symmetry. 



the scaling with size of deviation from particle-hole sym- 
metry and double-peaked nature of the order parameter 
histogram might suggest that the weak first order jump is 
tied intimately to the slight deviation from particle-hole 
symmetry. 

In order to test this hypothesis, we have located a se- 
ries of critical points along the phase boundary in the 
vicinity of the tip, and obtained detailed statistics of the 
charge-density wave order parameter for this sequence of 
critical points. If this hypothesis were valid, one would 
expect that the double-peak in the order parameter his- 
togram at the largest size would become successively less 
pronounced and extrapolate to zero upon approaching 
the tip of the lobe. Simultaneously, the deviation from 
quarter filling would also extrapolate to zero. 

However, our numerical results do not support this hy- 
pothesis. This is clear from Fig [131 where we show data 
at three critical points for our largest system size L — 48: 
Of these, one is hole-rich, the other particle rich, and the 
third is perfectly particle-hole symmetric to within our 
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FIG. 13: (Color online). Distribution of the order parameter 
m 2 and histogram of deviation of the density from quarter 
filling for L = 48, 0/L = 15/28. The legend indicates the 
value of / for each of these critical points. 



error bars. As is clear from the accompanying histogram 
of the charge-density wave order parameter, there is no 
significant smearing of the double-peak structure upon 
approaching the particle-hole symmetric transition. 

Thus, our numerical results lead us to conclude that 
although the transition at the tip appears to be 'close' 
to a landau forbidden second order transition, it is actu- 
ally a very weakly first order transition as far as we can 
tell with the largest sizes and lowest temperatures avail- 
able to us. As we now demonstrate below, an analytical 
study of the problem leads us to the conclusion that the 
effective low energy theory for the particle-hole symmet- 
ric transition is precisely the easy-plane NC CP 1 model of 
Refi^. Our conclusion of a weakly first-order transition 
then suggests that more work is needed to fully under- 
stand the physics of this universality class. 



IV. MAPPING TO A GAUGE THEORY 

We now develop a theory for the phase transition of the 
model in terms of the defects of the Mott insulator (MI) 





FIG. 14: (Color online). a)The four degenerate ground states 
associated with the plaquette VBS state, which may be asso- 
ciated with four different orientations of a Z4 order parameter 
ip c . Also shown in b) and c) are two types of Z4 vortices in the 
plaquette state. The blue dotted lines represent the four do- 
main walls. The encircled site is the core of the vortex. 



state at quarter filling, being careful to distinguish be- 
tween transitions with full low energy particle-hole sym- 
metry, and those without. The insulator phase is char- 
acterized by a Z4 clock order parameter ip c . In Fig [TJ] 
(a), we show the four plaquette ground states which may 
be associated with the four orientations of the Z4 order 
parameter (similarly, the four equivalent columnar states 
correspond to the four 'columnar' choices of the phase of 



Because the order parameter is discrete, the natural 
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topological defects are domain walls between different 
orientations of the phase of ip c , and four of them can 
terminate at a point to create a point-like Z4 vortex ex- 
citation. There are two types of vortices present, de- 
pending on whether the core of the vortex is hole-like 
with a negative charge deviation at that site (Fig [14] (b)) 
or particle-like with a positive charge deviation (Fig [T4l 
(c)). Translating the entire valence bond pattern by one 
lattice spacing reverses the direction of winding of the 
phase of the VBS order parameter. Thus, particle (hole)- 
like vortices have opposite windings on the two sublat- 
tices. Moreover, a particle-rich vortex has a phase wind- 
ing which is opposite to a hole-rich vortex located on the 
same sublattice. 

Keeping these two points in mind, we define the charge 
of a vortex with its core at site r*j of the underlying square 
lattice as n ri = (nu + U2% + n^i + na — 1) for r, G A 
sublattice, and n Ti = — (ni, + n2i + n^i + nn — 1) for 
Ti G B sublattice. Note that n Ti is non-zero only at the 
core of the vortices in any vortex configuration. Further, 
since the total boson charge may be expressed as 



Q = r ^ (nu + rin + n 3l + no) 



(5) 



we note that the vortices represent fractionally charged 
excitations carrying quanta of charge equal to one half 
the elementary boson charge. The mechanism for this 
is essentially the same as in many other examples dis- 
cussed earlier in the literature {e.g. the Motrunich- 
Senthil model realization of a Z2 deconfined phased). 
However, in contrast to earlier work on Z2 deconfined 
phase, these fractionally charged vortex excitations are 
expected to be strongly interacting in our model (see be- 
low), and we do not expect them to have an independent 
existence in the insulating phase. Nevertheless, as we 
will see below, they are convenient degrees of freedom 
in terms of which one may hope to describe the phase 
transition to a superfluid state, following the approach 
of Ref. [Tl| 

In addition, we define h TiT . as the charge on the links 
of the square lattice where , Tj are nearest neighbours 
on the underlying square lattice, and identify it with 
the boson number at the center of the link. We also 
find it convenient to soften the hard-core constraint on 
the boson number, and instead capture the same physics 
through a strong on-site Hubbard potential. Finally, it is 
also useful to work with a rotor description of the charge 
degrees of freedom (measured from a background value 
corresponding to the quarter filled insulating state). We 
can then describe the microscopic Hamiltonian in terms 
of rotor 'angular momenta' [n rA ,n rB ,h rArB ] , and their 
conjugate angle variables [ip rA ,ip rB1 (p rArB \. Writing the 
effective Hamiltonian in terms of these variables, we ob- 



tain 
H 



= Uj2(n rA + 1 - /) 2 + Uj2(-n rB + 1 - ff 

A B 

- 2t ^ C0S((^. A - lfr A , + ip rBrA ~ <fr B r A ,) 

r A r A , 

- 2t ^ C0S (¥V B ~ Vr B , + ¥r A r B , ~ <Pr A r B ) 

r B r B i 

+ VY,nl r , (6) 



along with the local constraints (due to the overcomplete 
set of variables used) 

^ t n rAr B i ~ 1 = n r A 



/ j > l r B r A < 



— 1 = — n r 



(7) 



As mentioned earlier, a large positive V is introduced 
to make the physics qualitatively similar to the hard-core 
boson case. This Hamiltonian and the local constraints 
define a compact U{1) gauge theory coupled to matter. 
To make this mapping explicit, we define a compact vec- 
tor potential a on the links of the square lattice such 
that a rArB = (f>r A r B and a TATB — —a TBTA . The conjugate 
electric field vector is then defined as e TATB = h r ATD , and 



tatb 1 

The Hamiltonian in these new variables 



reads 
H 



= Uj2(n rA + 1 - ff + Uj2(-n rB + 1 - ff 



A' 



2t cos((p rA , — ip rA — I a ■ dl) 

J A 
f B ' 

2t cos((p rB , — ip rB — / a - dl) 

r B r B , 



while the constraint can now be written as (V • e) r — 
n r + e r , where £a — 1 and cb = — 1- Thus, our effective 
model can be reformulated as a compact (a is an angle) 
17(1) lattice gauge theory coupled to two types of matter 
fields (n TA , n rB ) in the presence of a background (static) 
staggered charge pattern. 

The transformation properties of the fields (f>, n, a a , e a 
under discrete lattice symmetries are readily obtained 
(see Table lH|) and play an important role in the discussion 
below. The global ?7(1) symmetry associated with boson 
number conservation is realized by the phase rotation 

e** -> e iae "e^ (9) 

with a a constant. Correspondingly the total conserved 
boson number is given by 



Q 



(10) 
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Thus e xtpr A creates a state with total boson number 1/2 
while e llfr B creates a state with total boson number 
-1/2. 



Operation 


<t> 


n 


a a 


6q 


T x 


-<t> 


—n 


—a a 




Rtt/2 


-<t> 


—n 


—t-apap 




Ry 


-4> 


—n 


-a a 




T 


-d> 


n 


—a a 





TABLE II: Transformations of the fields on the square lattice 
under the various independent discrete symmetry operations. 
T x is translation by a unit cell in the i-direction. R w /2 is n/2 
rotation about the center of a plaquette. R y is reflection about 
the y-axis through midpoint of a horizontal bond. T is time 
reversal and is antiunitary. All other symmetry operations 
may be built up from these. 

This compact gauge theory is somewhat non-standard 
because the A type particles can never convert to the 
B type particles and vice-versa, which makes it a bit 
inconvenient to analyse. To overcome this technical dif- 
ficulty, we define a theory with two kinds of matter fields 
(rti,n2) (and their conjugate variables ((^1,(^2)), both of 
which can live on the A and B sublattices, and hop be- 
tween the A and B sublattice. The original restriction 
that particles of one type only live on the correspond- 
ing sublattice can be captured by a term of the kind 
(VJ2r B n i + ^Sr A n i) W1 "th V large and positive as be- 
fore. We expect that these changes preserve the universal 
properties of our original theory. 

Now, the charge at each site is (n\ T — n2r), and the 
constraint becomes (V • e) r = n\ r + n 2 r + e r . The elec- 
tric field is however still related to the microscopic boson 
density so that the electric field cor- 

relators give direct information about the density corre- 
lations in our microscopic model. Further simplification 
arises when we approach the transition at the tip of the 
MI lobe. Because of particle-hole symmetry at the tip, 
linear terms in ri\ , n 2 vanish. As a result it is enough for 
our purposes to consider the following simplified theory: 



- * / 5^cos(A itV »2-o M ) + V'5^^ 
- K J2cos(Ax a) + Vj2nl + Vj2nl 



□ r B 

(V • e) r = n>i r + n 2r + ( r 



(11) 



where n stands for the two lattice directions x and y. 
The K term has been added to make the theory look 
like a more standard gauge theory and represents a ring- 
exchange term in the dimer language. Such a term of this 
form would not change any of the universal properties of 



the theory, and is expected to be dynamically generated 
at low energies. The transformation properties of these 
new fields under lattice symmetries are also readily ob- 
tained and we list them below in Table As before 
global U(l) symmetry is realized by the phase rotation 



(12) 



Operation 


<f>3 


n s 


a a 


6q 


T x 


—4>s 


— Us 


—a a 




Rtt/2 


— <f>s 


— Us 






Ry 






—a a 




T 


-CPs 


n s 


—a a 





TABLE III: Transformations of the fields ((j> s ,n s ,a a ,e a ) on 
the square lattice under the various independent discrete sym- 
metry operations. Here s = 1,2 and s = 2 if s = 1 and 
vice- versa. 

This model may now be analysed by a duality transfor- 
mation. The results are rather similar to dual versions 
of boson models at half-filling on a square lattice. In 
particular there are two vortex species and ip2 (cor- 
responding to vortices in cpi and 4>2 respectively). These 
two vortices move independently on the dual square lat- 
tice. The universal critical properties of the transition 
may then be described in a continuum description of this 
dual theory. The oscillating background gauge charge 
(due to the e r in the Gauss law constraint) ensures that 
the leading term that allows for mixing of tpi and tp2 
vortices in the continuum dual action is 



d 2 xdr\ (^tp 2 ) 



c.c 



(13) 



The full action then reads 

s v 



J d 2 xdrJ2\( d ^- iA ^)^ 



+^l| 2 |</> 2 | 2 +u(|^l| 2 + |V'2| 2 ) 

+k (e Q/ 3 7 5^A 7 ) 2 



2 



(14) 



Here A^ is the dual gauge field whose magnetic field 
strength is the original boson density. This is exactly 
the same critical theory as that derived in Ref. [37] to 
describe the superfluid-VBS transition of bosons at half- 
filling on a square lattice. 

We thus expect that the present model lies in the same 
universality class as that transition. As argued in Ref. 
Hfi this is the easy plane NC CP 1 universality class and is 
described by the action S v . The term S m ix has been ar- 
gued to be irrelevant at the corresponding fixed point. 
This low energy theory S v was studied in Ref [38| by 
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Motrunich and Vishwanath. Their analysis concluded 
that the theory was self-dual and has a single continuous 
phase transition with unusual exponents that were esti- 
mated in their work based on a detailed numerical study 
of a particular lattice regularization of the theory. How- 
ever, in more recent work, Kuklov et. al. have studied 39 
a different lattice regularization and argued in favour of a 
weakly first order transition. Since the present work also 
finds a similar very weakly first order transition, more 
work appears to be needed to understand the physics of 
this universality class in greater detail. 

Finally, moving away from the tip of the Mott lobe 
corresponds to turning on a non-zero chemical poten- 
tial in the original boson model. The transition is then 
expected to be either first order or to proceed through 
two independent transitions with an intermediate super- 
solid phase. Although the present analysis cannot decide 
between these two possibilities without additional input 
specific to the microscopic model being considered, our 
numerics provides evidence in favour of a first order tran- 
sition. 



quantum phase transition from the superfluid state of the 
weakly interacting system to the lattice symmetry broken 
Mott insulator. Our analytical results above suggest a 
direct transition at the tip of the p — 1/4 insulating 
lobe — more precisely, this transition would be expected 
to be a continuous transition in the universality class of 
the easy-plane NCCP 1 model considered by Motrunich 
and Vishwanath in Ref l38l 

In our numerical work, we have studied the correspond- 
ing transition close to the tip of the lobe. Our numerical 
evidence points to a direct very weakly first order transi- 
tion close to the tip of the lobe. A detailed extrapolation 
of the strength of the first-order jump suggests that the 
transition remains very weakly first order even at the 
particle-hole symmetric tip. Thus, the present study un- 
derlines the need for further work to fully understand 
the physics of the easy-plane NCCP 1 universality class, 
particularly with a view towards understanding the very 
weakly first order nature of the transition seen in our 
work and in previous studies^ of models that belong to 
this universality class. 



CONCLUSIONS 
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We have thus obtained convincing numerical evidence 
for a four-fold lattice symmetry breaking Mott insulator 
at large boson repulsion in a simple model of hard-core 
bosons on the planar pyrochlore (checkerboard) lattice 
at quarter filling. We have also studied the nature of the 
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